Spatial Data as Matrices and Rasters

HES 505 Fall 2022: Session 10

Matt Williamson

Today’s Plan

Objectives

  • By the end of today, you should be able to:

    • Describe the raster data model and its representation in R

    • Access the elements that define a raster

    • Build rasters from scratch using matrix operations and terra

Defining the raster data model

The Raster Data Model

Defining the Raster Data Model

  • Vector data describe the “exact” locations of features on a landscape (including a Cartesian landscape)

  • Raster data represent spatially continuous phenomena (NA is possible)

  • Depict the alignment of data on a regular lattice (often a square)

    • Operations mimic those for matrix objects in R
  • Geometry is implicit; the spatial extent and number of rows and columns define the cell size

library(terra)
terra 1.7.39
filename <- system.file("ex/elev.tif", package="terra")
r <- rast(filename)
plot(r, col=gray(0:1000/1000), main = "Elevation (m)")

Types of Raster Data

library(stars)
Loading required package: abind
Loading required package: sf
Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
x = 1:5
y = 1:4
d = st_dimensions(x = x, y = y, .raster = c("x", "y"))
m = matrix(runif(20),5,4)
r1 = st_as_stars(r = m, dimensions = d)

r = attr(d, "raster")
r$affine = c(0.2, -0.2)
attr(d, "raster") = r
r2 = st_as_stars(r = m, dimensions = d)

r = attr(d, "raster")
r$affine = c(0.1, -0.3)
attr(d, "raster") = r
r3 = st_as_stars(r = m, dimensions = d)

x = c(1, 2, 3.5, 5, 6)
y = c(1, 1.5, 3, 3.5)
d = st_dimensions(x = x, y = y, .raster = c("x", "y"))
r4 = st_as_stars(r = m, dimensions = d)

grd = st_make_grid(cellsize = c(10,10), offset = c(-130,10), n = c(8,5), crs = st_crs(4326))
r5 = st_transform(grd, "+proj=laea +lon_0=-70 +lat_0=35")

par(mfrow = c(2,3), mar = c(0.1, 1, 1.1, 1))
r1 = st_make_grid(cellsize = c(1,1), n = c(5,4), offset = c(0,0))
plot(r1, main = "regular")
plot(st_geometry(st_as_sf(r2)), main = "rotated")
plot(st_geometry(st_as_sf(r3)), main = "sheared")
plot(st_geometry(st_as_sf(r4, as_points = FALSE)), main = "rectilinear")
plot(st_geometry((r5)), main = "curvilinear")
par(mfrow = c(1,1), mar= c(5.1, 4.1, 4.1, 2.1))

  • Regular: constant cell size; axes aligned with Easting and Northing

  • Rotated: constant cell size; axes not aligned with Easting and Northing

  • Sheared: constant cell size; axes not parallel

  • Rectilinear: cell size varies along a dimension

  • Curvilinear: cell size and orientation dependent on the other dimension

Types of Raster Data

  • Continuous: numeric data representing a measurement (e.g., elevation, precipitation)

  • Categorical: integer data representing factors (e.g., land use, land cover)

Continuous Rasters

mintemp <- rast("ftp://ftp.hafro.is/pub/data/rasters/Iceland_minbtemp.tif")
mintemp
class       : SpatRaster 
dimensions  : 340, 375, 1  (nrow, ncol, nlyr)
resolution  : 2000, 2000  (x, y)
extent      : -1050226, -300225.9, -699984.3, -19984.32  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=69 +lon_0=-4 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
source      : Iceland_minbtemp.tif 
name        : Iceland_minbtemp 
min value   :       -0.9982879 
max value   :        8.6031137 
plot(mintemp)

Categorical Rasters

  • Requires a classification matrix and coercion to factor

  • levels allows you to define the categories

# Create classification matrix
cm <- matrix(c(
  -2, 2, 0,
  2, 4, 1,
  4, 10, 2), ncol = 3, byrow = TRUE)

# Create a raster with integers
temp_reclass <- classify(mintemp, cm)
tempcats <- c("cold", "mild", "warm")
levels(temp_reclass) <- tempcats
Warning: [set.cats] setting categories like this is deprecated; use a two-column
data.frame instead
cats(temp_reclass)
[[1]]
  value category
1     0     cold
2     1     mild
3     2     warm
plot(temp_reclass)

Adding Dimensions

  • When data are aligned in space and/or time, more efficient to process as ‘cubes’ or ‘stacks’

  • Bands of satellite imagery, multiple predictors, spatio-temporal data

# (C) 2021, Jonathan Bahlmann, CC-BY-SA
# https://github.com/Open-EO/openeo.org/tree/master/documentation/1.0/datacubes/.scripts
# based on work by Edzer Pebesma, 2019, here: https://gist.github.com/edzer/5f1b0faa3e93073784e01d5a4bb60eca

# plotting runs via a dummy stars object with x, y dimensions (no bands)
# to not be overly dependent on an input image, time steps and bands
# are displayed by replacing the matrix contained in the dummy stars object
# every time something is plotted

# packages, read input ----
set.seed(1331)
suppressPackageStartupMessages(library(colorspace))
suppressPackageStartupMessages(library(scales))

# make color palettes ----
blues <- sequential_hcl(n = 20, h1 = 211, c1 = 80, l1 = 40, l2 = 100, p1 = 2)
greens <- sequential_hcl(n = 20, h1 = 134, c1 = 80, l1 = 40, l2 = 100, p1 = 2)
reds <- sequential_hcl(n = 20, h1 = 360, c1 = 80, l1 = 40, l2 = 100, p1 = 2)
purples <- sequential_hcl(n = 20, h1 = 299, c1 = 80, l1 = 40, l2 = 100, p1 = 2)
greys <- sequential_hcl(n = 20, h1 = 0, c1 = 0, l1 = 40, l2 = 100, p1 = 2)

# matrices from raster ----
# make input matrices from an actual raster image
input <- read_stars("img/slide_10/iceland_delta_cutout_2.tiff") # this raster needs approx 6x7 format
# if the input raster is changed, every image where a pixel value is written as text needs to be checked and corrected accordingly
input <- input[,,,1:4]
warped <- st_warp(input, crs = st_crs(input), cellsize = 200) # warp to approx. 6x7 pixel

# these are only needed for resampling
warped_highres <- st_warp(warped, crs = st_crs(warped), cellsize = 100) # with different input, cellsize must be adapted
# this is a bit of a trick, because 3:4 is different format than 6:7
# when downsampling, the raster of origin isn't so important anyway
warped_lowres <- st_warp(warped_highres[,1:11,,], crs = st_crs(warped), cellsize = 390)
# plot(warped_lowres)
# image(warped[,,,1], text_values = TRUE)

t1 <- floor(matrix(runif(42, -30, 150), ncol = 7)) # create timesteps 2 and 3 randomly
t2 <- floor(matrix(runif(42, -250, 50), ncol = 7))

# create dummy stars object ----
make_dummy_stars <- function(x, y, d1, d2, aff) {
  m = warped_highres[[1]][1:x,1:y,1] # underlying raster doesn't matter because it's just dummy construct
  dim(m) = c(x = x, y = y) # named dim
  dummy = st_as_stars(m)
  attr(dummy, "dimensions")[[1]]$delta = d1
  attr(dummy, "dimensions")[[2]]$delta = d2
  attr(attr(dummy, "dimensions"), "raster")$affine = c(aff, 0.0)
  return(dummy)
}

s <- make_dummy_stars(6, 7, 2.5, -.5714286, -1.14) # mainly used, perspective
f <- make_dummy_stars(6, 7, 1, 1, 0) # flat
highres <- make_dummy_stars(12, 14, 1.25, -.2857143, -.57) # for resampling
lowres <- make_dummy_stars(3, 4, 5, -1, -2) # for resampling

# matrices from image ----
make_matrix <- function(image, band, n = 42, ncol = 7, t = 0) {
  # this is based on an input image with >= 4 input bands
  # n is meant to cut off NAs, ncol is y, t is random matrix for time difference
  return(matrix(image[,,,band][[1]][1:n], ncol = ncol) - t)
  # before: b3 <- matrix(warped[,,,1][[1]][1:42], ncol = 7) - t2
}

# now use function: 
b1 <- make_matrix(warped, 1)
b2 <- make_matrix(warped, 1, t = t1)
b3 <- make_matrix(warped, 1, t = t2)
g1 <- make_matrix(warped, 2)
g2 <- make_matrix(warped, 2, t = t1)
g3 <- make_matrix(warped, 2, t = t2)
r1 <- make_matrix(warped, 3)
r2 <- make_matrix(warped, 3, t = t1)
r3 <- make_matrix(warped, 3, t = t2)
n1 <- make_matrix(warped, 4)
n2 <- make_matrix(warped, 4, t = t1)
n3 <- make_matrix(warped, 4, t = t2)

# plot functions ----
plt = function(x, yoffset = 0, add, li = TRUE, pal, print_geom = TRUE, border = .75, breaks = "equal") {
  # pal is color palette
  attr(x, "dimensions")[[2]]$offset = attr(x, "dimensions")[[2]]$offset + yoffset 
  l = st_as_sf(x, as_points = FALSE)
  if (li)
    pal = lighten(pal, 0.2) # + rnorm(1, 0, 0.1))
  if (! add)
    plot(l, axes = FALSE, breaks = breaks, pal = pal, reset = FALSE, border = grey(border), key.pos = NULL, main = NULL, xlab = "time")
  else
    plot(l, axes = TRUE, breaks = breaks, pal = pal, add = TRUE, border = grey(border))
  u = st_union(l)
  # print(u)
  if(print_geom) {
    plot(st_geometry(u), add = TRUE, col = NA, border = 'black', lwd = 2.5)
  } else {
    # not print geometry
  }
}

pl_stack = function(s, x, y, add = TRUE, nrM, imgY = 7, inner = 1) {
  # nrM is the timestep {1, 2, 3}, cause this function
  # prints all 4 bands at once
  attr(s, "dimensions")[[1]]$offset = x
  attr(s, "dimensions")[[2]]$offset = y
  # m = r[[1]][y + 1:nrow,x + 1:ncol,1]
  m <- eval(parse(text=paste0("n", nrM)))
  s[[1]] = m[,c(imgY:1)] # turn around to have same orientation as flat plot
  plt(s, 0, TRUE,  pal = purples)
  m <- eval(parse(text=paste0("r", nrM)))
  s[[1]] = m[,c(imgY:1)]
  plt(s, 1*inner, TRUE,  pal = reds)
  m <- eval(parse(text=paste0("g", nrM)))
  s[[1]] = m[,c(imgY:1)]
  plt(s, 2*inner, TRUE,  pal = greens)
  m <- eval(parse(text=paste0("b", nrM)))
  s[[1]] = m[,c(imgY:1)]
  plt(s, 3*inner, TRUE, pal = blues) # li FALSE deleted
}

# flat plot function
# prints any dummy stars with any single matrix to position
pl = function(s, x, y, add = TRUE, randomize = FALSE, pal, m, print_geom = TRUE, border = .75, breaks = "equal") {
  # m is matrix to replace image with
  # m <- t(m)
  attr(s, "dimensions")[[1]]$offset = x
  attr(s, "dimensions")[[2]]$offset = y
  # print(m)
  s[[1]] = m
  plt(s, 0, add = TRUE, pal = pal, print_geom = print_geom, border = border, breaks = breaks)
  #plot(s, text_values = TRUE)
}

print_segments <- function(x, y, seg, by = 1, lwd = 4, col = "black") {
  seg = seg * by
  seg[,1] <- seg[,1] + x
  seg[,3] <- seg[,3] + x
  seg[,2] <- seg[,2] + y
  seg[,4] <- seg[,4] + y
  segments(seg[,1], seg[,2], seg[,3], seg[,4], lwd = lwd, col = col)
}

# time series ----

# from: cube1_ts_6x7_bigger.png
offset = 26
plot.new()
#par(mar = c(3, 2, 7, 2))
par(mar = c(0, 0, 0, 0))
#plot.window(xlim = c(10, 50), ylim = c(-3, 10), asp = 1)
plot.window(xlim = c(-15, 75), ylim = c(-3, 10), asp = 1)
pl_stack(s, 0, 0, nrM = 3)
pl_stack(s, offset, 0, nrM = 2)
pl_stack(s, 2 * offset, 0, nrM = 1)
# po <- matrix(c(0,-8,7,0,15,3.5,  0,1,1,5,5,14), ncol = 2)
heads <- matrix(c(3.5, 3.5 + offset, 3.5 + 2*offset, 14,14,14), ncol = 2)
points(heads, pch = 16) # 4 or 16
segments(c(-8, 7, 0, 15), c(-1,-1,3,3), 3.5, 14) # first stack pyramid
segments(c(-8, 7, 0, 15) + offset, c(-1,-1,3,3), 3.5 + offset, 14) # second stack pyramid
segments(c(-8, 7, 0, 15) + 2*offset, c(-1,-1,3,3), 3.5 + 2*offset, 14) # third stack pyramid
arrows(-13, 14, 72, 14, angle = 20, lwd = 2)  # timeline
text(7.5, 3.8, "x", col = "black")
text(-10, -2.5, "bands", srt = 90, col = "black")
text(-4.5, 1.8, "y", srt = 27.5, col = "black")
y = 15.8
text(69, y, "time", col = "black")
text(3.5, y, "2020-10-01", col = "black")
text(3.5 + offset, y, "2020-10-13", col = "black")
text(3.5 + 2*offset, y, "2020-10-25", col = "black")

A note about support

  • We talked briefly about the agr option in the st_sf() function

  • agr refers to the attribute-geometry-relationship which can be:

    • constant = applies to every point in the geometry (lines and polygons are just lots of points)
    • identity = a value unique to a geometry
    • aggregate = a single value that integrates data across the geometry
  • Support is the area to which an attribute applies.

  • Rasters can have point (attribute refers to the cell center) or cell (attribute refers to an area similar to the pixel) support

Rasters in R

Rasters in R

  • raster - the original workhorse package; built on sp, rgdal, and rgeos
    • RasterLayer, RasterStack, and RasterBrick classes
  • terra - relatively new; developed by the raster folks, but designed to be much faster
    • SpatRaster and SpatVector classes
  • stars - developed by sf package developers; tidyverse compatible; designed for spatio-temporal data
    • stars class
    • Crosswalk between raster and stars is available here
    • Only way to deal with rectilinear and curvilinear data

Rasters with terra

  • syntax is different for terra compared to sf

  • Representation in Environment is also different

  • Can break pipes, Be Explicit

Rasters by Construction

mtx <- matrix(1:16, nrow=4)
mtx
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16
rstr <- terra::rast(mtx)
rstr
class       : SpatRaster 
dimensions  : 4, 4, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 4, 0, 4  (xmin, xmax, ymin, ymax)
coord. ref. :  
source(s)   : memory
name        : lyr.1 
min value   :     1 
max value   :    16 
plot(rstr)

Note: you must have raster or terra loaded for plot() to work on Rast* objects; otherwise you get Error in as.double(y) : cannot coerce type 'S4' to vector of type 'double'

Rasters by Construction: Origin

  • Origin defines the location of the intersection of the x and y axes
r <- rast(xmin=-4, xmax = 9.5, ncols=10)
r[] <- runif(ncell(r))
origin(r)
[1] 0.05 0.00
r2 <- r
origin(r2) <- c(2,2) 
par(mfrow = c(2,1))
plot(r, xlim = c(-1, 6), ylim=c(0,7), main = "Original")
plot(r2, xlim = c(-1, 6), ylim=c(0,7), main = "New Origin")

par(mfrow=c(1,1))

Rasters by Construction: Resolution

  • Geometry is implicit; the spatial extent and number of rows and columns define the cell size
  • Resolution (res) defines the length and width of an individual pixel
r <- rast(xmin=-4, xmax = 9.5, 
          ncols=10)
res(r)
[1] 1.35 1.00
r2 <- rast(xmin=-4, xmax = 5, 
           ncols=10)
res(r2)
[1] 0.9 1.0
r <- rast(xmin=-4, xmax = 9.5, 
          res=c(0.5,0.5))
ncol(r)
[1] 27
r2 <- rast(xmin=-4, xmax = 9.5, 
           res=c(5,5))
ncol(r2)
[1] 3

Rasters from Files

  • Building rasters useful for templates
  • More common to read from files
r <- rast(system.file("ex/elev.tif", package="terra"))
r
class       : SpatRaster 
dimensions  : 90, 95, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : elev.tif 
name        : elevation 
min value   :       141 
max value   :       547 

Accessing Raster Attributes: Coordinate Reference System

  • terra stores CRS in WKT format

  • Can set and access using EPSG and proj (deprecated)

  • Pay attention to case

r <- rast(system.file("ex/elev.tif", package="terra"))
crs(r)
[1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"

Accessing Raster Attributes: Coordinate Reference System

r <- rast(system.file("ex/elev.tif", package="terra"))
crs(r, describe=TRUE)
    name authority code area         extent
1 WGS 84      EPSG 4326 <NA> NA, NA, NA, NA
crs(r, proj=TRUE)
[1] "+proj=longlat +datum=WGS84 +no_defs"
crs(r, parse=TRUE)
 [1] "GEOGCRS[\"WGS 84\","                                   
 [2] "    DATUM[\"World Geodetic System 1984\","             
 [3] "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,"   
 [4] "            LENGTHUNIT[\"metre\",1]]],"                
 [5] "    PRIMEM[\"Greenwich\",0,"                           
 [6] "        ANGLEUNIT[\"degree\",0.0174532925199433]],"    
 [7] "    CS[ellipsoidal,2],"                                
 [8] "        AXIS[\"geodetic latitude (Lat)\",north,"       
 [9] "            ORDER[1],"                                 
[10] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"
[11] "        AXIS[\"geodetic longitude (Lon)\",east,"       
[12] "            ORDER[2],"                                 
[13] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"
[14] "    ID[\"EPSG\",4326]]"                                

Accessing Raster Attributes: Bounding box

  • terra uses ext() to get or set the extent/bounding box

  • Fills cells with NA

ext(r)
SpatExtent : 5.74166666666667, 6.53333333333333, 49.4416666666667, 50.1916666666667 (xmin, xmax, ymin, ymax)
r2 <- r
ext(r2) <- c(5, 7, 48, 52)
ext(r2)
SpatExtent : 5, 7, 48, 52 (xmin, xmax, ymin, ymax)
plot(ext(r2))
plot(ext(r), add=TRUE, col="blue", main="Blue is original extent")

Converting vectors to rasters

  • Sometimes necessary to convert between data models

  • raster::rasterize, terra::rasterize, stars::st_rasterize, and fasterize::fasterize all will convert polygons to raster data

  • stars::st_polygonize will work in the opposite direction

  • terra::vect will read in vectors as SpatVectors or coerce sf to SpatVector